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A variety of swimming microorganisms, called ciliates, exploit the bending of a large number 

of small and densely-packed organelles, termed cilia, in order to propel themselves in a viscous 

fluid. We consider a spherical envelope model for such ciliary locomotion where the dynamics of the 

individual cilia are replaced by that of a continuous overlaying surface allowed to deform tangentially 

to itself. Employing a variational approach, we determine numerically the time-periodic deformation 

of such surface which leads to low-Reynolds locomotion with minimum rate of energy dissipation 

(maximum efficiency). Employing both Lagrangian and Eulerian points of views, we show that in the 

optimal swimming stroke, individual cilia display weak asymmetric beating, but that a significant 

symmetry-breaking occurs at the organism level, with the whole surface deforming in a wave-like 

fashion reminiscent of metachronal waves of biological cilia. This wave motion is analyzed using a 

3 formal modal decomposition, is found to occur in the same direction as the swimming direction, and 

^ is interpreted as due to a spatial distribution of phase-differences in the kinematics of individual 

cilia. Using additional constrained optimizations, as well a constructed analytical ansatz, we derive 

a complete optimization diagram where all swimming efficiencies, swimming speeds, and amplitude 

' ' of surface deformation can be reached, with the mathematically optimal swimmer, of efficiency one 

I I half, being a singular limit. Biologically, our work suggests therefore that metachronal waves may 

f5 allow cilia to propel cells forward while reducing the energy dissipated in the surrounding fluid. 

I 

5^ I. INTRODUCTION 

. ^ Swimming microorganisms are found in a large variety of environments. From spermatozoa cells to bacteria in 

^ the human body, from unicellular protozoa to multicellular algae swimming in the ocean, these small organisms are 
able to generate net locomotion by exploiting their interaction with a surrounding viscous fluid Jj . Because of their 
small dimensions, the Reynolds number associated with their motion, Re = UL/v, where U is the swimming velocity, 
L the typical size of the organism, and v the fluid kinematic viscosity, is close to zero and the effects of body and 
fluid inertia are both negligible. As a result, the motion of swimming microorganisms is based on the exploitation of 
viscous drag to generate thrust through periodic non-time-reversible shape changes [THS|. 

Many swimming microorganisms use the beating of elongated flexible appendages attached to their surface to 
^\ produce motion 0]. These appendages are known as cilia or flagella depending on their distribution density on 
7—1 the cell, and their size relative to that of the organism. Eukaryotic flagella, such as those used by invertebrate 

^NJ or mammalian spermatozoa, are longer than the cell head and are only found in small numbers (typically one for 

I — [ spermatozoa, and a few for other eukaryotes [T]). In contrast, most of the surface of ciliates such as Paramecium (see 
f~^ Fig. Ill left) is covered by cilia much shorter than the cell body [5] . Ciliary motion is also functionally essential for 
f^ respiratory systems, where the beating of cilia covering lung epithelium permits the transport of mucus and foreign 
^-H particles out of the respiratory tract [6 . In that case, the cilia support is fixed and the cilia motion produces a net 

flow. 

Eukaryotic flagella and cilia have similar diameters (of the order of 200 nm), beating frequencies (of the order of 10 
Hz) and internal structure. They differ however in length. The typical length of a sperm cell flagellum is of the order 
of 50 ^m, while cilia, such as those of Paramecium, have a typical length of 5-10 //m W. Both eukaryotic flagella 
and cilia are subject to distributed actuation through the sliding of neighboring polymeric filaments (microtubules 
doublets) [II H]. Notably, bacterial flagella, although bearing the same name, are much smaller (typically 20 nm 
diameter and 5 — 10 /xm in length), have a much simpler internal structure than their eukaryotic equivalent, and are 
passively driven by a rotating motor located in the cell wall [HI HHI ■ 

The propagation of wave patterns is essential to both flagellar and ciliary propulsion [1 Ij . The deformation of 
individual eukaryotic flagella in planar or helical wave patterns leads to non-time-reversible kinematics, and is thus 
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FIG. 1: Paramecium swimming using metachronal waves (Hoffman modulation contrast microscopy, 40X magnification). Left: 
Picture of tlie cell in a microcfiannel. Riglit: Visualization of metachronal waves of cilia deformation propagating along the 
cell; time increases from top to bottom, with a time difference of 1/300 s between each picture; arrows indicate the propagation 
of a wave of effective strokes in the cilia array. Scale bar is 50 fim in both pictures. Pictures courtesy of Sunghwan (Sunny) 
Jung, Virginia Tech. 



responsible for the net swimming motion of the associated microorganism [T]-[31 [TH [13] • In this paper, we focus on 
cihary propulsion, for which two levels of symmetry-breaking (and two types of waves) are observed |^ . At the level 
of an individual cihum, a deformation wave propagates along the cilium length as for flagellar motion. The beating 
pattern is however different from that of individual flagella, and the individual stroke of a cilium can be decomposed 
into two parts: an effective stroke, during which the cilium is extended and offers the most resistance to the fluid, 
and a recovery stroke, in which the cilium is bent in such a way as to reduce the viscous drag. 

In addition to such asymmetric beating at the level of an individual cilium, the beating coordination of neighboring 
cilia at the organism level results in a collective behavior known as metachronal waves. All cilia on the surface of a 
microorganism perform similar beating patterns, but they deform in time with a small phase difference with respect 
to their neighbors, and these phase differences are spatially-distributed in a way that leads to symmetry-breaking at 
the level of the whole cell, and the formation of a wave pattern of surface deformation IHHlj (see Fig. [l] right) . The 
origin of the synchronization responsible for the metachronal waves in ciliary propulsion is still debated, but several 
recent studies have suggested that it results from hydrodynamic interactions between neighboring cilia |15H19| . 

In this paper, we consider the energy cost and hydrodynamic efficiency associated with ciliary propulsion. By 
efficiency, we understand here a relative measure of the organism displacement or velocity to the energy dissipated 
through viscous stresses in the flow to produce this motion. Of course, this differs from the actual energetic cost for 
the organism which includes metabolism and other internal biological considerations, and we limit ourselves here to 
the purely hydrodynamical aspect of the efficiency. Although little is known experimentally about the actual energy 
consumption associated with ciliary propulsion, some studies suggest that metachronal waves reduce the energy loss 



PP] , More generally, the question of the hydrodynamic efficiency of different swimming patterns is at the heart of 
many low-Re locomotion investigations of the extent to which the swimming modes observed in nature could be 
optimal with respect to hydrodynamic efficiency. In particular, numerous studies have considered the optimal beating 
of flagella [2J I21H24| . Through a theoretical and numerical optimization framework, we determine in this paper the 
particular collective cilia beating patterns that minimize the dissipation of mechanical energy. 

Computing the flow around a ciliated organism accurately is difficult because of the large number of appendages 
deforming and interacting. Two types of modeling approaches have been proposed in the past to address this problem. 
The first type of model, termed sublayer modeling, considers the dynamics of individual cilia, either theoretically in 
a simplified fashion [U [S] , or numerically with all hydrodynamic interactions [33 [5S] . The study of the collective 
dynamics of a large number of cilia remains however a costly computation jl5[ I20j . 

An alternative approach, motivated by the densely packed arrangement of the cilia on the surface of the organism, 
is based on the description of the swimmer by a deformable, continuous surface enveloping the cilia at each instant 
[H [271 [21] ■ This so-called envelope model substitutes the motion of material surface points (cilia tips) with the 
deformation of the continuous surface, and is expected to be a good approximation when the density of cilia is 
sufficiently high. In the particular case of a purely spherical swimmer, this model is known as a squirmer |281I29| . and it 
was recently used to study the collective dynamics and the rheology of suspensions of model swimming microorganisms 
[30l[3T] . Although some organisms using cilia do have a spherical shape (e.g. Volvox [32] ) , most ciliated microorganisms 
have an elongated body. Nonetheless, the squirmer model and the spherical approximation allow one to reduce the 
complexity of the problem in order to shed some light on the fundamental properties of symmetry-breaking in ciliary 
locomotion. 

In this paper, we thus consider locomotion by a squirmer in a Newtonian fluid without inertia as a model for 
locomotion of a spherical ciliated cell. We will assume that the squirmer can deform its shape tangentially in a 
time-periodic fashion, and hence the shape remains that of a sphere for all times. For a given time-periodic stroke — 
that is for a given periodic Lagrangian surface displacement field — we are able to compute the swimming velocity of 
the swimmer, the energy dissipation in the fluid, and use both to define the stroke hydrodynamic efficiency. We will 
not restrict our analysis to small-amplitude deformations 1331 134| , but will allow arbitrary large-amplitude tangential 
deformations to take place. The purpose of our work is to then determine the optimal squirming stroke which 
maximizes this swimming efficiency (or, alternatively, minimizes the amount of work done against the fluid for a fixed 
average swimming velocity), and to study its characteristics. In a first part, we will consider the general problem 
of determining the optimal stroke theoretically and numerically. In a second part, we will consider the limitation 
to the surface kinematics introduced by the finite size of the cilia, and how it affects the optimal swimming stroke. 
In all cases, we will show that the optimal swimming strokes show very little asymmetry at the level of individual 
cilia (Lagrangian framework) but display strong symmetry-breaking at the level of the whole organism (Eulerian 
framework) , reminiscent of phase differences between cilia and of metachronal waves observed in biology. Within our 
framework, we will thus be able to conclude that metachronal waves are hydrodynamically optimal. 

In Sec. [nj after a brief review of the squirmer model and its dynamics, the optimization procedure is presented, 
togetherw with the resulting optimal stroke. In Sec. |III[ a constraint is added to the optimization to take into account 
the finite-length of the cilia and limit the surface displacement. Based on the observations of the optimal strokes 



obtained in Sees. [TI] and III an analytical ansatz is constructed in Sec. jIVj that achieves asymptotically the theoretical 
upper bound for the swimming efficiency of a squirmer. In Sec. [Vl the physical properties of the optimal strokes are 
presented and discussed, in particular the wave characteristics. Our results are finally summarized and discussed in 
Sec.[Vll 

II. OPTIMAL SWIMMING STROKE OF SPHERICAL SWIMMER 

We consider in this paper the dynamics of a spherical microorganism able to produce locomotion by imposing time- 
periodic tangential displacements of its spherical surface — the so-called squirmer approximation — as an envelope 
model for ciliated cells. Only axisymmetric strokes with no azimuthal displacements are considered, thereby restricting 
the swimming motion to a pure time-varying translation. 

A. Equations of motion and swimming efficiency 

1. Geometric description 

The dynamics of the micro-swimmer is studied in a translating reference frame centered at the squirmer center. In 
this co-moving frame, the fluid and surface motions are described in spherical polar coordinates (r, 6, ip) (see notation 
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FIG. 2: Notation for the squirmer model. A co-moving frame centered on the squirmer center is used with spherical polar 
coordinates {r,9,if). The axisymmetry of the surface motion imposes a purely translating swimming motion along e^ with 
time- varying velocity U{t). See text for details. 



in Fig. I2|. All quantities are non-dimcnsionalized using the swimmer radius a, the stroke frequency / and the fluid 
dynamic viscosity fi. 

The surface of the squirmer (r — 1) is described by two Lagrangian variables < ^o !i ^r and < ipo l£ Stt, and 
the position of each surface element is defined by the time evolution of its polar and azimuthal angles 0{Oq, ipo, t) and 
ip(9Q,ipQ,t). Considering purely axisymmetric deformations of the swimmer's surface, we have 



0^i9{eo,t), ifi^cpo. 



(1) 



In the following, we will therefore omit the dependence in ipQ, all particles located on a circle of constant 9 having 
the same behavior. By a straightforward symmetry argument, the swimming velocity of the organism must thus be 
of the form U{t)ex. 

Each point on the swimmer surface can be equivalently described by its polar angle 9 or its vertical cartesian 
coordinate x = cos 9. In the following, both notations will be used to describe the surface motion 



= t9(6lo,t), or X ^£,{xQ,t). 



(2) 



Most of the analytical work will be done in terms of the x variable as equations adopt a simpler form in this case. 
Physical properties of the stroke will however be discussed using the variable 9 which describes the actual angular 
displacement of the cilia tips on the surface. 

The Lagrangian label xq = cos^o is chosen such that a;o = — 1 (^o == ^) and xq = 1 (^o — 0) are respectively the 
south and north poles of the squirmer (Fig. ^. Spatial boundary conditions at both poles, as well as time-periodicity 
impose that 



e(-l,t) = -l and e(l,i) = 
^{xo,t) = £.{xo,t + '2Tr) for 



1 for < i < 27r, 
- 1 < Xo < 1. 



(3) 

(4) 



Note that there exist multiple choices for the Lagrangian label xq (or 9q), and it is not necessarily the position of 
the material point at t — 0. Any other Lagrangian label xi satisfying Eq. ([3]) can be used provided that Xq — >■ xi is 
a bijective function. In the following, we use the mean value of d over a swimming period as Lagrangian label and 
^0 = (^)i where (•) denotes the averaging operator in time. 

Finally, in a discrete representation of the three-dimensional motion of cilia, it would possible for cilia to cross. This 
is however forbidden in the continuous axisymmetric envelope model used here, and the surface velocity is uniquely 



defined at each point; ^(a;o) niust ttierefore be a bijective function, which leads to 



dxo 



> 0, for - 1 < a;o < 1, < t < 2-k. 



(5) 



In the co-moving frame, the surface velocity u'^ is purely tangential u"^ — Ug{6,t)e0, and is related to the surface 
displacement by a partial time derivative 






{df,,t)^ul{^{9o,t),t). 



Equivalently, this equation can be rewritten using the variables {x, t) as 



dt 



{xo,t) — u{^{xo,t),t), with u{x,t) = —\/l — x'^Ug {cos ^ x,t). 



(6) 



(7) 



The tangential velocity Ug{9,t) is the velocity at a fixed point 6, while d'd/dt is the tangential velocity of a given 
material point labelled by 6*0. Equation ^ is therefore the fundamental conversion from Eulerian to Lagrangian 
quantities. Note that u{x,t) in Eq. ^ is the axial component of the surface velocity (along e^). 



2. Swimming m,otion of a squirmer 

In the limit of zero Reynolds number Re ~ Q, the equations for the incompressible flow around the squirmer simplify 
into the non-dimensional Stokes equations 



V^u = Vp, V.u = 0, 

with boundary conditions, expressed in the co-moving frame, as 

u = u'^ = u^ {6, t)eg forr = l, 
u = —UGx for 7' — > CX3. 



Eqs. ([8|)-(10l can be solved explicitely as ;2S 
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where Ln(x) is the 77-th Legendre polynomial, Vn{9) is defined as 



(2r7 + l)sin0 

L„(cost'), 



77(77+ 1) 



(8) 
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(11) 
(12) 
(13) 

(14) 



and the time coefficients a„ {t) are obtained uniquely from the expansion of the surface tangential velocity in spherical 
harmonics 



ue(e,t)=uf(0,t)=E""(*)^"(^)- 



(15) 



In the Stokes regime, the inertia of the swimmer is negligible so the fiuid force is zero at all times 

Text = / [-pn + (Vu + Vu"^) • n] d5 = 0. (16) 

Js 



The viscous torque is also zero by symmetry. From Eqs. (11)-(13| and (16), the swimming velocity U as well as the 
instantaneous rate of work of the swimmer on the fluid can then be obtained as 1281 



U{t)^ai{t), 



V{t) = 127r 
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(17) 

(18) 



The swimming velocity is completely determined by the first mode ai{t), which is thereafter referred to as the 
swimming mode. All others modes, a„(i), do not contribute to the swimming velocity but do contribute to the energy 
consumption, and in that sense are penalizing the swimming efficiency of the organism. However, the existence of 
these non-swimming modes is required to ensure the periodicity of the surface displacement. 



3. Swimming efficiency 



In this paper we are going to derive the optimal stroke kinematics, and thus we have to define our cost function. 
Since low-i?e locomotion is essentially a geometrical problem |13l[35], the appropriate cost function is basically a way 
to normalize this geometrical problem. Here we will use the traditional definition of a low-_Re swimming efficiency, ry, 
given by 
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(19) 



where I'* is the force required to drag a rigid body of same shape as the swimmer at the time-averaged swimming 
velocity (t/), and {V) the mean rate of energy dissipation in the fluid during swimming [TH51 [^I5S] . The cost function, 
77, given in Eq. (19), has traditionally been termed an efficiency, but perhaps more accurately it should be termed a 



normalization, as it is the ratio between the average rate of work necessary to move the body in two different ways: 
in the numerator, dragging the body with an external force, and in the denominator, self-propelled motion at the 
same speed. It is thus not a thermodynamic efficiency, and for general swimmers does not have to be less then one, 
although for most biological cells it is on the order of 1% [T]. 

In the particular case of a squirmer, Stokes' formula gives T* = 67r([/), and the efficiency (Eq. 19) can be obtained 
analytically as 
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(20) 



From the Cauchy-Schwartz inequality, we see that Eq. (20) leads to an upper bound of 77 < 1/2. This upper bound. 



which we will discuss in more detail in Sec. IV is tighter than the one obtained in Ref. |34j (77 < 3/4) and corresponds 
to the treadmilling microswimmer described in Ref. [37 . Importantly, the definition of efficiency retained here is 
purely mechanical and does not characterize the absolute efficiency of the locomotion mode for the organism, which is 
influenced by many other factors, including feeding and metabolic costs. In addition, since we are using an envelope 
model for the cilia, our approach does not allow us to capture the fluid dissipation in the sublayer (i.e. near the cilia). 
Using Eqs. ^ and (15), we have the dynamics 
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Multiplying both sides of the previous equation by L'^ (^) and integrating in ^ leads to 
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Integrating by part and using the boundary condition at the poles (Eqs. pi flnally leads to the Lagrangian-Eulerian 
relationship 
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(23) 



For a given stroke x — ^{xo,t), the coefficients a„(i) are computed from Eq. (23) and used to determine the 
swimming efficiency 77 [^] using Eq. (20). 



B. Optimization of the swimming efficiency 

The objective of the present work is to determine the stroke ^(xo,i) maximizing the swimming efficiency 77[f]. 
Physically, the optimal stroke will therefore be the one swimming the furthest for a given amount of dissipated 
energy, or, alternatively, will be the one minimizing the rate of work done by the swimmer for a given swimming 
speed. To characterize the optimal stroke, we use a variational approach. 



1. Swimming efficiency gradient in the stroke functional space 

Let us consider a particular reference swimming stroke ^{xo,t), and the perturbed stroke ^{xo,t) + d^{xQ,t), where 
SS.{xo,t) is a small perturbation to the reference stroke. In the following, tilde quantities correspond to the reference 
(or, initial) stroke. From the boundary conditions (Eq. pi) applied to ^ and ^ + 5^, we obtain the following boundary 
conditions and periodicity constraints on 6S,{xQ,t) 



S^{-1, t) = 6^(1, t) = for all < t < 27r, 
S^{xo,t) ^S^{xo,t + 2Tr) for ah - 1 < xo < 1. 

Retaining only the linear contribution, the change in swimming efficiency is given by 
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From Eqs. (23) and (24 1, the perturbation (5a„(i) induced on a„(t) is computed as 
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For n = 1 and using Eq. ( 25 ) , the time average of Eq. ( 28 1 leads to 
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Keeping only the leading order contribution, we have 



Using Eq. (28) and integration by part in time, we find 
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After substitution into Eq. ( 26 ) , the variation Sr/ can finally be written as 

/I p2TT 
J F[^]{x,m{xo,t)dxodt, 

where F[^](xo, t) is the efficiency gradient in the stroke functional space evaluated at f 
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(24) 
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(29) 
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(32) 



(33) 



2. Projection on the subspace of the boundary conditions 
F[£] is the gradient of ij with respect to the swimming stroke f . Choosing a sufhciently small S£, aligned with F[£] 



guarantees an increase of the swimming efficiency {6ri > 0) between S. and ^ + (5^ as the integrand in Eq. (32) is positive 



everywhere, which suggests the manner in which we can numerically iterate to find the optimal swimmer. However, 
such a change S^ along F[^] does not guarantee that C + <5? will satisfy the boundary and periodicity conditions (Eqs. 
p}|4| or the constraint on monotoneous variations (Eq. pi) . To circumvent this difficulty, we notice that for any stroke 
£,{xo,t), there exists at least one continuous function ilr(xo,t) such that we can write 



^(xo.t) 
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J^^mx',tfdx' ' 



(34) 



Performing the optimization on the field il'ixo, t) rather than (,{xq, t) frees us from imposing the boundary, periodicity 
and bijection conditions separately, as, through Eq. (34), any 4'i^O: t) will generate an acceptable £,{xo, t). The analysis 



of the previous paragraph remains valid and we now must evaluate the variation 5£,{xo,t) induced by a small change 
S4>{xo,t). From Eq. (34), we have 
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Substitution of this result into Eq. (32) leads, after integration by part and rearrangement of the result, to 

r2TT 
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/o J-i 
with G[ip]{xo,t) the efficiency gradient in the functional space of "ip, that can be written as: 
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with F[£,] given in Eq. (33) 



C. Numerical optimization 



The optimization on the swimming stroke is performed iteratively and numerically, using a steepest ascent algorithm 
|38) . Given a guess ^nixo,t), the efficiency r][tpn] and the gradient G[V'n] are computed numerically using Eqs. (20l, 
(|23[), (33) and (37). The next guess ■0„_|_i(a;o, t) is obtained numerically by marching in the functional space in the 



gradient direction as 



V'„+i(a;o,t) = ipnixo,t) + eG[ipn]{xo,t), 



(38) 



where £ is a small number. For small enough e, Eq. (36) guarantees that Sr] — rjlijjn+i] ~ vH'n] > 0. Starting from an 
initial guess i/'Oi s-nd using this approach iteratively, one travels in the (ip^rj) -space along the steepest slope in rj until 
convergence is reached. As the purpose of the present paper is to obtain physical insight on the maximum efficiency 
strokes, we did not attempt to construct the fastest converging algorithm, and it is possible that faster convergence 
would be achieved, for example, using a variation of the conjugate gradient algorithm. 

The efficiency 77 and gradient G are computed numerically from ip using spectral methods in both t (Fast Fourier 
Transform) and xq (Chebyshev spectral methods). To avoid aliasing phenomena introduced by the successive non- 
linear products, the 2/3-dealiasing rule is applied in both space and time before taking each physical product. 

Convergence is achieved when marching in the direction of the gradient G does not induce any increase of the 
swimming efficiency anymore, even after successive reductions of the step size e. This algorithm, as most iterative 
optimization techniques, can not guarantee the finding of an absolute maximum of the swimming efficiency but only of 




FIG. 3: The contribution of the successive modes, (a^), to the efficiency (stars) decreases exponentially with the mode order, 
n. The results are plotted for the unconstrained optimization with efficiency rj ~ 22.2%, and were obtained with M — 40, 
N^ = f50 and Nt = 128. 



local maxima. Because of the infinite number of dimensions of the functional space, it is expected and observed that 
some optimization runs will lead to local maxima with small values of ?/. This difficulty is overcome by performing 
several runs with different initial conditions and/or resolution until convergence to similar solutions provide enough 
confidence in the finding of an absolute maximum. 

As pointed out before, S,{xq, t) is not a unique description of the swimming stroke as any bijection on the lagrangian 
label will lead to another equivalent representation of same 77. To identify whether two strokes are equivalent, we 
must therefore resort to the comparison of their physical characteristics or velocity distribution. 

We use Nx Chebyshev Gauss-Lobatto points in xq, and Nt equidistant points in t. The infinite sum in Eq. (20) 



must be truncated at M modes. A large enough number of modes M must be computed for the efficiency estimate 
to be accurate, thereby imposing a practical minimum on N^ and Nt for the integrals involved in Eq. ( 23 ) to be 



computed accurately. In the representative cases presented in the rest of this work, typical values of the discretization 
parameters were M « 40-80, N^ « 120-200 and Nt « 64-128. Figure^ shows that the contribution (a^) of mode n 
to the efficiency of the reference stroke decreases exponentially with n. One can then estimate the truncation error 
on 77 as « 0.03%. 



D. Initial guesses for the swimming stroke 



Different sets of initial conditions were used to ensure convergence of the procedure, 
conditions are 



Examples of such initial 



^{xQ,t) = xo + ci{l-Xg) cos{t-C2Xo), 



(39) 



with N an integer, ci ^ 1 and 0.5 < C2/7r < 10. Eq. (39) corresponds to a traveling wave in xq of wavespeed l/c2, 
attenuated at both poles. The efficiency being a quadratic function of the swimming velocity (see Eq. 20), the initial 
condition must have non-zero swimming velocity and must therefore not be time-reversible (Purcell's theorem |13j). In 
the following, traveling wave patterns will be identified and it is important to ensure that the wave-like characteristics 
of the optimal solution result from the efficiency optimization and not from the particular initial guess considered. In 
that regard, different wave amplitudes and wave velocities were tested in Eq. (39). Also, a superposition of multiple 
waves traveling in different directions were tested as well as "hemispheric solutions" defined independently in each 
hemisphere (the equator then becoming a fixed point), and "elliptical solutions" of the form 



^{xo,t) 



(l + a;o) 



1 + ci sin ( ^^^ 



(40) 
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FIG. 4: Lagrangian description of the optimal swimming stroke with no constraint on the maximum displacement. The 
maximum amplitude of displacement is 0max ~ 52.6°, the mean swimming velocity is (f/) ~ 0.33 and the swimming efficiency 
is »7 « 22.2%. Each curve illustrates the trajectory 9 = ^{9o,t) of a single material surface point. 



Independently of the particular choice of initial conditions, we obtain a systematic convergence to the results presented 
in the following section. 



E. Unconstrained results: optimal swimming stroke 



For all these different initial conditions, we observe the convergence of our numerical approach toward a stroke of 
efficiency i] w 22.2%. We will refer to this solution as the optimal unconstrained stroke thereafter. 

The trajectories of surface elements along the sphere 9 = i?(0o,O ^i'^ shown on Fig. El and Fig. [5 presents a 
sequence of snapshots corresponding to a stroke period. On Fig. [5] the color code is a Lagrangian label and allows one 
to track the position of a particular surface particle in time. A few lines Oq — constant are also represented for better 
visualization of the surface motions. The average swimming velocity associated with this stroke is directed upward 
and equal to (U) = 0.33. Due to the length and time scaling chosen here, this corresponds to a mean swimming 
velocity slightly above one body-length-per-period. The swimming velocity is not constant throughout the period as 
can be seen on Fig. [6J 

The optimal unconstrained swimming stroke illustrated in Figs. |4] and [5] can be decomposed, in time, into two 
parts: (1) an effective stroke where the surface moves downward (increasing 9) while stretching; (2) a recovery stroke 
where the surface elements that migrated toward the south pole during the effective swimming stroke are brought 
back toward the north pole, in the same direction as the swimming velocity. The surface is highly compressed in this 
phase, which gives it a shock-like structure (see the dark region in Fig. B|. The instantaneous swimming velocity, 
U{t), is maximum and roughly constant throughout the effective stroke, while the recovery stroke is associated with 
a reduced (even reversed) swimming velocity (Fig. pi). 

These two strokes are not exactly successive in time as their boundary is not vertical on Fig. [4] The recovery stroke 
is mostly visible during the middle half of the stroke period (7r/2 < t < 3n/2) (see Figs. Illandls]). Outside of this 
time domain, the entire organism surface is moving downward (effective stroke), with the possible exception of two 
small regions located near the poles. This is confirmed by the predominance of the swimming mode ai{t) over the 
other modes outside of the domain [tt/2 , 37r/2] (see Fig.l6|. 

Although individual cilia motion is not explicitly represented in this continuous envelope model, the swimming 
stroke obtained through this optimization process shares many similarities with metachronal waves observed in ciliated 
microorganisms 01 JHI ■ Indeed, a small phase difference in the motion of neighboring surface points (or cilia tips) can 
be observed: as can be seen in Fig.|4J for the range 20° < 9q < 160°, the maximum and minimum of a given trajectory 
occur with a small time delay compared to those of a trajectory with a slightly greater value of ^o- 
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FIG. 5: (color online) Snapshots of the displacement of the squirmer surface over a swimming period in the optimal uncon- 
strained stroke (Omax ~ 52.6°, (U) « 0.33, ry « 22.2%). The surface color of each material point on the surface refers to the 
mean polar angle of that material point. Black lines also correspond to the location of particular material points, and have 
been added for clarity. 



This small phase difference results in a global wave pattern at the organism level propagating from the south pole 
{9 — 180°) to the north pole {9 — 0°) as clearly seen on Fig. [5] We illustrate schematically in Fig. [Tlthis symmetry- 
breaking mechanism through the collective behavior of individual Lagrangian points. With purely identical and 
symmetric motions of equal amplitude of neighboring Lagrangian points, the introduction of a small spatial phase 
shift breaks the symmetry and generates a shock-like dynamics similar to the recovery stroke obtained numerically 
and shown in Fig. [4] For biological cells, the biophysical origin of the observed phase-locking between neighboring cilia 
and the generation of metachronal waves is still a matter of investigation f l5H19j . Our results show that a wave- like 
deformation of the surface at the whole-organism level is actually an optimal for the swimming efficiency. 

Finally, we note that the efficiency obtained through our optimization procedure is much larger than that of typical 
swimming microorganisms. Our optimal unconstrained stroke is found to be associated with very large displacements 
of material points on the spherical surface, with a maximum amplitude of ©max ~ 52.6° (meaning that there exists 
a Lagrangian point for which 'd{9o,t) — 9q varies between —52.6° and -1-52.6°), corresponding to a linear tangential 
displacement equivalent to 90% of a diameter. This stroke of large amplitude is however not realistic for real ciliated 
organisms as the length of the cilia (and therefore the maximum distance covered by its tip) is typically several 
times smaller than the size of the organism itself ^3^ . In the following section our optimization approach is adapted 
to constrain the maximum displacement amplitude of each surface point, and we study the maximum achievable 
swimming efficiency with a given maximum individual tangential displacement. The optimal unconstrained stroke 



12 




FIG. 6; Time-evolution of the ten first modes, a„(t) {n < 10), over the period for the optimal unconstrained stroke shown in 
Fig. [5] The swimming mode or swimming velocity U{t) — ai{t) is plotted as a thick solid line. Modes 2 through 10 are plotted 
as dashed line. The swimming mode is observed to dominate during the effective stroke (for t < 7r/2 and t > 37r/2), while 
many modes are significant during the recovery stroke (7r/2 < t < 3n/2). 





(a)In-phase beating 



(b)Small phase difference 



FIG. 7: A small phase-difference between trajectories of neighboring Lagrangian points with identical and symmetric individual 
trajectories can lead to symmetry breaking and the appearance of a shock-like structure, (a) Harmonic in-phase beating, (b) 
Same harmonic motion with a linearly spatially varying phase. 



obtained above will serve as a reference for comparison and discussion in the rest of the paper. 



III. STROKE OPTIMIZATION WITH CONSTRAINED SURFACE DISPLACEMENTS 
A. Applying a constraint on the maximum displacement 



The optimization algorithm described in Sec. |TT]can be adapted to include a constraint on the maximum displace- 
ment of individual points on the surface, by modifying the function to be maximized as 



j^m-j H 



{T[^\{xo) ~c)Axq, 



(41) 



13 



180 




FIG. 8: (color online) Same as Fig. 14] in a constrained-displacement case. Here, the maximum displacement is 0max ~ 10.6° 
the mean swimming velocity is (U) ~ 0.058 and the swimming efficiency is ?7 « 6.2%. 



where rj [^] is the efficiency as defined in Eq. ( 20 1 , T [^] is a measure of the amplitude of the displacement of individual 
material points for the stroke ^{xQ,t)^ and c is a dimensionless threshold parameter (a smaller c corresponding to a 
stricter constraint). H is defined as 



H{u) = A 



1 



tanh I — 



u\ 



(42) 



This form of H, when A is large and e is small, introduces a numerical penalization in the cost function only when 
the displacement measure T is greater than the threshold value. Values of A = 10"* and e = 10""^ were typically used 
in our numerical calculations, with little effect on the final result. 
Physically, the constraint T[f] = Q[£] should ideally be applied, with 



mK^o) 



.(Oo) 



'^t [cos 1 {^{xQ,t))\ -mint [cos ^ {^{xo,t))\ 



(43) 



which is the actual displacement amplitude of an individual point. However, the strong non-linearity of this measure 
is not appropriate for the computation of a gradient in functional space as presented in Sec. Illj An alternative measure 
of the displacement is 



T[e]M = ((e-(e»^ 



1 



S,{xQ,t) 



1 

2^ 



axo,t')dt' 



di, 



(44) 



which is the variance of the displacement along the vertical axis. Note that in using Eq. (44), we do not strictly 



enforce that the maximum displacement along the surface should be less than a given threshold. In addition, choosing 
Eq. (44) instead of Eq. (43) introduces a difference between points located near the poles (a small change in 9 there 
corresponds to a much smaller vertical displacement compared to that obtained with the same angular displacement 
at the equator). This is however not an issue as Eq. (44) will still penalize large amplitude strokes and the constraint 
applied using Eq. (44) will be stronger near the equator, where the displacements were observed the largest in the 
unconstrained problem (Fig. HI). In addition, a posteriori, an implicit relationship between c and the actual maximum 
displacement can be obtained, and c can thus be used as a tuning parameter to constrain the maximum amplitude 

Following the same approach as in Sec. Hj we consider a small perturbation 5^ of a reference swimming stroke ^. 
The resulting change in the cost function J[^] is obtained as 



6J ^ 



2tx 



^[i]{xo, t)6£,{xo, t)dxodt, 



(45) 
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FIG. 9: Same as Fig. [5] in the constrained case. The example illustrated here has a maximum displacement amplitude 
6max ~ 10.6°, a mean swimming velocity (U) « 0.058 and a swimming efficiency r] « 6.2%. 



where the modified gradient J^ is now 



m = m--{c-{o]H' T[^]- 



(46) 



and the corrected gradient G[ip] of the cost function with respect to ■0 is then obtained from T as in Eq. (37). 



B. Results 



A number of simulations were performed using a range of values of c, each with various choices for both A and 
e. As expected, we observe that introducing the penahzation and reducing the value of c does indeed reduce the 
maximum displacement amplitude of individual Lagrangian points. This decrease in 0max is found to systematically 
be associated with a decrease in the swimming efficiency. The swimming stroke obtained under constrained maximum 
displacements presents the same structure as the unconstrained case (see Figs, [s] and |9]). The {0,t) domain can be 
divided into two strokes, an effective stroke where the surface moves downward while stretching, and a recovery stroke, 
corresponding to the upward motion of a compressed surface. The maximum amplitude of displacement is however 
smaller, as is the propagation velocity of the shock structure corresponding to the recovery stroke (Fig. Isl). Because 
of this smaller shock velocity, the recovery and effective strokes co-exist over the entire swimming period on different 
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Case 


(a) (b) (c) (d) (e) (f) (g) 


fe'max 


52.6° 45.1° 34.2° 26.9° 14.4° 11.0° 4.26° 


vi%) 


22.2 21.0 17.9 15.1 8.33 6.81 2.61 


(U) 


0.33 0.29 0.22 0.17 0.075 0.051 0.02 



FIG. 10: Variations of the swimming velocity over a period for optimal strokes with different maximum displacements. The 
unconstrained optimal stroke is represented by a thick black line (a). The values of the maximum displacement, the swimming 
efficiency and the mean swimming velocity are also given for each case. 



regions of the swimmer's surface, resulting in a smoothing of the swimming velocity who becomes close to a constant 
when the constraint on Omax becomes tighter (Fig. 10). 

For each simulation, the maximum ampUtude 0rnax — niax(0[^](a;o)) can be computed together with the swimming 
efficiency rj. A clear monotonic relationship between rj and 0max is obtained, as shown in Fig. 11 (left). The dispersion 
around the main trend is small and due to some numerical runs converging to local minima. The right-most point 
of the curve corresponds to the unconstrained optimization presented in Sec. |ll] Similarly, the swimming velocity is 
found to be an increasing function of the maximum tangential displacement in the optimal stroke (Fig. 11 right). It 
is maximum for the unconstrained stroke (0max ~ 52.6°) with a dimensionless speed of [/ « 0.33, corresponding to 
about one body length per period, and is reduced to less than half a body length per period when the displacement 
amplitude is constrained below sa 30°. 



IV. THEORETICAL UPPER BOUND ON THE SWIMMING EFFICIENCY OF A SQUIRMER 

Using the numerical optimization approach described above, we have obtained an optimal swimming stroke with 
an efficiency of ~ 22.2%. We have further shown that constraining the maximum displacement of the surface reduces 
the swimming efficiency continuously from that maximum value. In this section, we are addressing the question of 
whether the numerical result of about 22.2% is the maximum achievable efficiency for a periodic swimming stroke. 
We show below that, in fact, the theoretical upper bound of 50% can be reached asymptotically by a singular stroke. 



A. Bound on the swimming efficiency and conditions of equality 



The Cauchy-Schwartz inequality states that (al) > (ai)^ and the equality can only be reached when ai{t) — (ai) 



is a constant. From the definition of the swimming efficiency ry in Eq. (20 1, we therefore obtain that for any swimming 
stroke £_ 



V[^]< 



1 



(47) 
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FIG. 11: (color online) Left: Swimming efficiency, 77, vs. the maximum amplitude in angular displacement of the stroke, Omax- 
The results of multiple runs (with different initial conditions or penalizing function parameters) have been regrouped into bins. 
Horizontal and vertical error bars represent the standard deviation in displacement and efficiency for each point. The blue 
star corresponds to the optimal unconstrained swimming stroke of Sec. Ill] Right: Swimming velocity vs. maximum angular 
displacement. 



and this upper bound can only be reached for the particular case: 

ai(t) — (ai) and a„(i) = for all n > 2. 



(48) 



In this case, ue(9,t) ^ sm6 and the vertical surface velocity is constant and uniform. This stroke corresponds to 
the so-called "treadmilling" swimmer |37) which is not time-periodic as the poles are permanent source and sink of 
material surface points, which move continuously from one pole to the other. 

However, the treadmill swimmer provides some important insight about strokes maximizing the efficiency: ai must 
be dominant over the other modes to minimize the energy consumption. This observation is consistent with the 
relative mode amplitudes an{t) observed in the result of the optimization procedure (Fig. p|. To guarantee the stroke 
periodicity, at least one of the two conditions in Eq. (48 1 must be violated. In the following, we build an analytical 



ansatz that achieves 77 = 1/2 asymptotically based on the observation that the conditions of Eq. (48) need only to be 
satisfied on most of the period, the interval where they are violated being asymptotically of zero measure. 



B. Building the analytical ansatz 

The analytical ansatz is the superposition of two parts, the effective stroke or outer solution which covers most of 
the (x, i)-plane, and the recovery stroke or inner solution, which takes the form of a shock and enforces periodicity of 
the surface displacement. 



1. Outer solution 



The outer solution satisfies the optimality conditions of Eq. ( 48 ) exactly. The Lagrangian equation of motion of 
the surface points is obtained from Eqs. ([t]) and (15 1 as 



— =7(x -1) with 7=-ai, 



(49) 



which can be integrated as 



^{xo, t) = tanh ( log \/ ^ _ '^° - -ft 



(50) 



17 



180p^-r^Tr-=; 




(a)Trajectories 9 = •d{9o,t) and shock position 
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(b)Trajectories z = (;{zo,t) and shock position 



FIG. 12: The analytical ansatz is decomposed into an outer solution (effective stroke) in which only the swimming mode is 
present and an inner solution (recovery stroke) of small width which enforces periodicity in time. The trajectories corresponding 
to the outer solution are shown in dashed lines, and thick solid lines show the position of the shocks and inner solution trajectories 
(Eq. 54 1. Here, 7=1 and /3 = 0.7, and three swimming periods are shown. 



Using the change of variables x = tanh^/ and ^ ~ tanhx, Eqs. (15) and (50l become 

dx f . v^ 2n + 1 / N ^/ r w \i 

-Kriyo^t) = - > —r — — 7vQn(^)4i tanh(x) , -00 < y < 00, -tt < t < tt, 

at ^-—^ n n + 1) 

Outer solution: x(yo5 1) = yo — 7^, —00 < yo < 00. 
The trajectories of the outer solution are straight lines in the (x, t)-plane and are not periodic since 

x(2/o, tt) - x(yo, -tt) = -277r. 



(51) 
(52) 

(53) 



2. Shocks and inner solution 
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To construct a time-periodic solution, thin shock-like structures are super imposed to this outer solution (Fig. 12 1 
We look for shock equations of the form: 



xs„ = tanh[/3(i + 2fc7r)] or ys^ = /3(f + 2/c7r). 



(54) 



In the (y, i)-plane, a Lagrangian surface particle will follow a straight line of slope —7, until it hits one of the shocks 
(Fig. 12). The velocity inside the shock must be designed in such a way that the Lagrangian point will reemerge 



from the shock when Ay = 2ttj. The following change of variables, z = y and t = t — y/ (3 stretches the horizontal 
coordinate and positions the shocks as vertical lines r = 2A;7r in the (z, T)-plane (see Fig. [12]d). The trajectories of 
Lagrangian particles y = xiuoj t) ^re now rewritten as z = (^(zq, t), with 



dx di, 



1 



(3dT 



dt dr 

In (z, r) coordinates, the outer solution trajectories can be expressed as 

7T 



<r(20,T) = Zq 



I+7//3 



(55) 



(56) 



The Eulerian velocity w{z,t) in (z, r)-coordinates is now a function of r only (Fig. 12) and 'w{t) = —^fiji^^ + /?) 
everywhere, except in thin regions near t = 2A;7r. Considering the Gaussian approximation for (5cr(T), 



the solutions of 



dT 



SAr) 



w{t) with w{t) 



(Jy/n 



7 + /3 



(27r Y^ 5^{T + 2k-K)-l) 



(57) 



(58) 



/c— — 00 



are periodic and equal to the outer solution outside thin regions of characteristic width a. Changing back to variables 
(?/,i) and (a;,i), the ansatz is finally obtained as 



dt 



u{x,t) = (l-a;2) 



E 'A'- 



-7 + (/? + 7)- 



fc — — CO 



tanh X 



2kTr 



27r7 



+ E ^A'- 



k— — oo 



tanh X 



/3 



2fc7r 



(59) 



For given values of /3, a and 7, the efficiency of the corresponding periodic stroke can be computed numerically 
by projecting u{x,t) using Eq. (22). The infinite sums in Eq. ([59]) can be easily truncated after a few terms as the 



contribution of larger values of k to the term in brackets is limited to the vicinity of ±1. 



C. Asymptotic convergence to rj = 50% 



The solution of Eq. ( 59 ) satisfies the periodicity constraint and matches the (optimal) treadmilling swimming stroke 
over most of the domain. For a small value of the shock width a, we numerically sweep the parameter space (/3,7), 
and display the map of the corresponding values of the swimming efficiency in Fig. |13[ W e observe that when (3 and 
7 become very large, values of the efficiency approaching 50% can be reached. Figure 13 also shows the variations of 
77 with CT for large values of f3 and 7. We thus obtain that the value of 77 = 50% can be obtained in the asymptotic 
limit where cr — > and /3, 7 — > 00. 

The first criterion {a — >■ 0) is expected to occur as the shock must be narrow for the measure of the domain where 
the outer solution does not hold to be small. In addition, we see that the shock velocity (/3) and swimming velocity 
(27/3) must also go to infinity. This is the result of the non-dimensionalization of the equations using the stroke 
frequency. When 7 — > 00, the effective stroke accumulates an asymptotically infinite amount of surface at one pole 
between two successive recovery strokes; /? — > 00 corresponds to the recovery stroke taking as little time as possible 
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FIG. 13: (color online) (Left) Variations of the swimming efficiency obtained for the analytical ansatz, for a sweep in the (/3, 7) 
parameter space, and a shock width parameter of a = 0.05; lines display iso-values of rj. (Right) Evolution with a of the 
efficiency of the analytical ansatz for 7 = 1000 (solid) and 7 = 200 (dashed); in both cases, /3 — 800. 



to bring back this accumulated surface near the opposite pole, thereby to a time-periodic stroke. If the swimming 
velocity had been used for non-dimensionalization, the limit 7 — > (X) would have corresponded to an infinite spacing 
between two successive recovery strokes. The present asymptotic solution with /?, 7 — ^ cx) therefore matches the 
treadmill swimmer except for a time-interval of measure zero. In this singular (and thus, non-physical) limit, the 
entire surface of the sphere travels from one pole to the other during the effective stroke before the shock/recovery 
stroke redistributes the surface points, corresponding to a maximum displacement amplitude Qnmx is equal to 90°. 



D. Efficiency vs. maximum displacement — Final optimization diagram 



We therefore established numerically that for any value of rj less than, but arbitrarily close to, 50%, one can find 
parameter values for /3, 7 and a for which the stroke (Eq. 59) leads to swimming with efficiency rj. This upper bound 



is therefore reachable asymptotically but the corresponding stroke is singular. It is therefore not surprising that such 
a stroke, where all Lagrangian points accumulate at the south pole, could not be obtained through the numerical 
optimization procedure presented in Sees. |TT]and|TlT] which is based on a Lagrangian description of the surface 

The analytical approach can however also be used to confirm the results of our optimization approach for small 
values of the swimming efficiency. Equation (59) defines a family of periodic strokes with three parameters /3, 7 and a. 



For each stroke, we can compute numerically its efficiency together with the maximum displacement amplitude 6max- 
The envelope curve, or equivalently the maximum efficiency obtained for a given maximum displacement amplitude. 



can then be compared to the results of the optimization procedure from Secs.Oand III The results of this comparison 
are displayed in Fig. [M] 

In the range 0° < Qamx l£ 52°, solutions can be obtained both with the analytical ansatz and with the numerical 
optimization. We see that the maximum efficiency obtained with the ansatz stroke agrees with the optimal curve 



obtained in Sec. Ill confirming the validity and form of the bounds on 77 introduced by a restriction on 0max in 
our numerical approach. Beyond the value 0max — 52°, the analytical ansatz provides an indication of the bound 
imposed on the efficiency by a weaker restriction on Omaxj in a domain not reachable by our numerical algorithm, 
and allows us to continuously link the results of the numerical optimization results to the theoretical optimal, but 
singular, stroke of 77 = 1/2. 

Figure [T4| presents thus the complete optimization diagram, where for each value of 6max, the stroke has been 
optimized to lead to the largest possible swimming efficiency, and all values of 77 from to 1/2 can be obtained. In 
particular, we observe that realistic values of the cilia-length-to-body-length ratio lead to the angular amplitude range 
©max ~ 5-15°, which corresponds to optimal strokes of efficiency in the range 77 « 3-6% comparable to, if not slightly 
above, swimming efficiencies expected using flagellar propulsion [1] |4] . 
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FIG. 14: (color online) Final optimization diagram: Optimal efBciency of the swimming stroke, t), as a function of the maximum 
Lagrangian angular displacement of the surface, 0max- Green squares: result from numerical optimization (the vertical and 
horizontal error bars denote the variability from the initial conditions and the weight of the constraint); blue star: unconstrained 
optimal swimming stroke with 0max ~ 52.6° and 77 ~ 22.2°; black dots: envelope of the results obtained from the analytical 
ansatz for a — 0.1 and varying /3 and 7; red dot: asymptotic limit of the singular quasi-treadmill swimmer of efficiency ry — 1/2 
and maximum angular amplitude 0max ~ 90°. 



PROPERTIES OF OPTIMAL STROKES 



In this section, we consider the optimal strokes obtained numerically using the unconstrained and constrained 
optimization algorithms of Sees. [TI|and|III[ and analyze their physical properties. 



A. Lagrangian vs. Eulerian quantities 



Two points of view can be adopted to analyze the surface deformation of an optimal squirmer. One can either 
consider the property of a fixed point in the swimmer's co-moving frame characterized by a polar angle 9 (Eulerian 
formulation) or the property of a given material surface point indexed by ^o (Lagrangian formulation) . A quantity Q 
can then be understood and plotted either as a Lagrangian quantity Q{0o, t) or an Eulerian quantity Q{0, t). To make 
this distinction, 9 will be used in what follows for the fixed coordinate, and ^o will be used as a Lagrangian label. 

The mean position of a given material point ^o = (^) is one of the many possible Lagrangian labels. It is physically 
convenient and intuitive because it is expected to correspond roughly to the position of the cilia base. Note however 
that the configuration 'd{9o,t) = 9q for all ^oi corresponding to the case where all cilia are vertical at the same time, 
is never reached during the stroke period. 

In our analysis, we will see that quite different conclusions can be reached when considering the same quantity 
from one point of view or the other. This is illustrated below by looking at the stretching of surface elements. The 
Lagrangian formulation allows to address the individual behavior of a given material point or cilia tip in time. In 
contrast, the Eulerian formulation is more adapted to study the global properties of the swimming stroke, determined 
by the collective behavior of the different surface points. 
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Optimal unconstrained 
n ^ 22.2%, e„,ax ~ 52.6°, {[/) = 0.33 




Optimal constrained 
n ~ 6.2%, Omax « 10.6° , (17) = 0.058 
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FIG. 15: (color online) (Top) Lagrangian and (bottom) Eulerian map in the (6, t) plane of the stretching parameter 5 for the 
unconstrained swimming stroke (left) and the result of constrained optimization with efficiency 6.2% and maximum displacement 
10.6° (right). 



B. Surface stretching and compression 



The difference in the displacement of two neighboring surface elements induces some periodic stretching and com- 
pression of the surface. This stretching can be quantified using S — dii/ddo, with 6*0 the mean position of each surface 
point; S > I (resp. S < 1) corresponds to a stretched (resp. compressed) surface in comparison with the reference 
(mean) configuration. The stretching field S can be plotted either in Lagrangian coordinates as S{9o, t) or in Eulerian 
coordinates as S{6, t), and both are shown in Fig. 15 for two different strokes (the optimal unconstrained stroke with 
77 ss 22.2% and the constrained one with ry w 6.2%). 

One can see from the comparison of the top and bottom figures in Fig. [15] that the Lagrangian or Eulerian point 
of view fundamentally changes the description of surface stretching. In the Lagrangian formulation, the stretching of 
a material particle appears symmetric, with a particle spending about as much time compressed than stretched. In 
the Eulerian approach, however, the conclusion is qualitatively different. In this case, a given location experiences 
a stretching of the surface for most of the period, with compression being only achieved during the fast passage 
of the recovery shock. We therefore observe a strong symmetry-breaking at the organism level in stretching and 
compression, whereas it is essentially symmetric at the individual cilium level. This symmetry-breaking, consistent 
with the discussion presented in Sec. |IIE[ will be analyzed quantitatively in Sec. |VC[ 

Note in addition that in Fig. |15[ by comparing the stretching maps for the unconstrained and constrained strokes. 
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we can observe additional changes such as a smaller velocity of the shock in the constrained case. These properties 
will be discussed further in Sec. IV Dl 



C. Symmetry-breaking and collective behavior 

1. Individual (Lagrangian) motion asymmetry 



We first adopt the Lagrangian point of view to characterize the symmetry or asymmetry in the motion of individual 
surface elements. Quantitatively, different measures of this asymmetry can be chosen and we focus here on the 
following ones: 

(1) temporal asymmetry, measured as the period fraction during which the material point is moving downward (resp. 
upward) 



A 



Ci(0o)=log 



(S.o.)>0)\ 



yr,-(^o,t)<0 



(60) 



(2) maximum velocity asymmetry, defined as the ratio between the maximum velocities of the material point in each 
direction 



C2(0o)=log 



Min,( 1^(00, t) 



f)f) 



(61) 



(3) individual trajectory eccentricity, defined as the distance of the turning points to the mean position 



C3(eo)=log 



70 — "mill 



(^o) 



(62) 



These three quantities are plotted as functions of the Lagrangian coordinate, ^o, on Fig. [16] A logarithmic scaling 
is used so that (^p = corresponds to the symmetric situation and opposite sign ^p corresponds to the same amplitude 
of asymmetry in opposite directions. Regardless of the chosen asymmetry measure, we see that it vanishes near the 
poles, and the equatorial points correspond to local minima of the asymmetry. The maximum individual asymmetry 
is reached around 9o ~ 40-45°, and this maximum is more pronounced for C2 and C3 measuring a displacement 
asymmetry than for the temporal asymmetry factor, d. We also see that C,i and C2 are positive everywhere (except 
near the poles, where the very small displacements make the asymmetry measurements irrelevant), while C3 has 
opposite sign depending on the hemisphere, showing an eccentricity of the individual trajectories directed toward the 
equator. 

The value of the temporal asymmetry, Ci, is rather small and uniform for individual particles (Ci ~ 0.2 corresponds 
to a 20% asymmetry in the velocity sign) confirming the near-symmetry observed in the map of surface velocity 
displayed in Fig. [TBk. Individual cilia spend thus about the same time moving in either direction, but a slight 
asymmetry is observed, with about 20% more time spent performing the effective stroke (positive 86 /dt). As a 
consequence, the maximum velocity achieved is greater during the recovery stroke (positive ^2)- 



2. Group (Eulerian) asymmetry 



Although individual material points display a near-symmetric motion over the whole period, the collective behavior 
of the different surface elements and the phase difference in their motion leads to a symmetry breaking at the 
organism level. This was illustrated above by considering the stretching/compression parameters S in both Lagrangian 
and Eulerian coordinates, and observing that the near-symmetric Lagrangian patterns show significant Eulerian 
asymmetry. 



A similar disparity is apparent between the Lagrangian (Figs. 16 1) and Eulerian (Figs. 17 d) tangential velocity 



maps. The Eulerian asymmetry can be quantified using the normalized temporal asymmetry ratio Ci, the Eulerian 
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FIG. 16: (color online) Asymmetry in the motion of individual surface points (Lagrangian approach), (a) Lagrangian velocity 
map, ug{6o,t). (b) Asymmetry in the motion of individual material points measured as: the ratio of the period fraction with 
upward or downward velocity i^i (solid), the ratio of the maximum (in absolute value) positive and negative velocities ^2 (dash- 
dotted) and the excentricity of individual trajectories (3 (dashed). All three measures are plotted for the unconstrained stroke 
(r? ^ 22.2% and Gmax ~ 52.6°). 
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FIG. 17: (color onl ine) (a) Surface velocity map of the Eulerian velocity, ug{6,t). (b) Eulerian asymmetry measured as the 
ratio Ci(^) (Eq. 63 1 between the period fractions where the velocity at a fixed point is directed upward vs. downward. Data 
plotted for the optimal unconstrained stroke (r] ~ 22.2% and Omax ~ 52.6°). 



equivalent to C,i (Eq. 60): 



CiW^log 



T{u{e,t) >o) 

T{u{6,t) <0) 



(63) 



and the results are shown on Fig. [TTJ j. In comparison to Fig. [TBj j, we observe a significant increase of this asymmetry 
factor in comparison to the individual asymmetry measurement. While individual material points spend only a 
maximum of 20% more time performing the effective vs. the recovery stroke, a fixed Eulerian location on the 
swimmer near the equator is part of the effective stroke for 90% of the period. Note that while the results of Figs. 16 



and 17 are shown for the unconstrained optimal stroke, similar behavior is observed for the strokes with constrained 
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displacement amplitudes. 

In optimal strokes, although the motion of individual points on the swimmer surface is roughly symmetric, the 
collective behavior of all surface points display a global symmetry-breaking between the two phases of the swimming 



stroke (effective and recovery). As discussed in Sec. IV such breaking of the front-back symmetry in the squirming 
dynamics is essential to efficient swimming in order to ensure that the effective stroke dominated by the swimming 
mode covers most of the {6, i)-plane. 

D. Modal decomposition of the optimal stroke and w^ave propagation 

The emergence of global asymmetry from near-symmetric individual motion typically occurs when the individual 
displacement is harmonic with a spatially- varying phase. For example, trajectories of the form 

9 = ^{00, t)=eo + a{eo) cos(i - Oo/ve) (64) 

would correspond to exactly symmetric individual trajectories, but for sufficiently large amplitude a(0o), non- 
symmetric collective behavior is observed, thereby breaking the symmetry between upward and downward strokes 
(see Fig. [7]). 

1. Complex Empirical Orthogonal Function (CEOF) decomposition 

The symmetry-breaking identified in the optimal strokes is associated with a propagating wave reminiscent of the 
metachronal waves observed on the surface of many ciliated organisms ^. Here, we confirm the existence of this 
wave and quantify its properties through the use of a complex modal decomposition [39]. This method is based on 
the extension of classical Empirical Orthogonal Functions (EOF - also known as Proper Orthogonal Decomposition) 
to complex variables, and is well-designed to study propagating two-dimensional wave patterns |40j . 

In the classical EOF method, a function of space and time, /(6'o, t), is decomposed in a series of modes 

n 

designed in such a way that truncation of the infinite sum at any order N minimizes the L2-iiorm of the remainder. 
The Complex EOF decomposition performs the classical EOF decomposition on the Hilbert transform F{9o,t) of 
f{OQ,t). This technique is used here to study the displacement of surface points from their mean position 9q, which 
is decomposed as 



d{0o,t)-9o =Re 



Y,MOo)Bn{t) 



J2 an{Oo)bn{t) cos(0„(f?o) - Mt))- (66) 



This method is particularly powerful for our time-periodic problem. For all the optimal swimming strokes obtained 
in Sees. O and III the first mode is so dominant that the error introduced by restricting the sum in Eq. (|66[) to the 



first mode is found to be less than 1-2%. The characteristics (temporal and spatial amplitude, and phase functions) 



of this first mode are shown in Fig. 18 for both the unconstrained and constrained optimal strokes with decreasing 



displacement amplitude. Note that the definition of a„ and 6„ in Eq. ( 66 1 is not unique, and one is free to choose a 
normalization for one of these functions. We chose to normalize the temporal amplitude so that (bi) — 1. Similarly, 
the phase functions (pi and tpi are defined up to an arbitrary additive constant. 

2. Propagating wave mode and phase velocity 



In all cases plotted on Fig. 18 the temporal phase ?Ai(i) is found to be a linear function of time with unit slope, 
which is consistent with the 27r-stroke period. The spatial phase 4>i{9q) varies also linearly, except in narrow regions 
near the poles. Outside these polar regions, the first CEOF mode is therefore a progressive wave of phase velocity 

ye = ^iM, (67) 
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FIG. 18: Amplitude and phase characteristics of the dominant Complex Empirical Orthogonal Function (CEOF) mode in the 
unconstrained swimming stroke of maximum efficiency (thick solid hne - 0max ~ 52.6° and -q ~ 22.2%) and in constrained 
optimal strokes (thin lines): Omax ~ 44° and 77 ~ 21% (dashed), 0max ~ 25° and rj ~ 12.5% (dash-dotted), 0max ~ 11° and 
77 « 6.8% (dotted). 



where dot and prime denote respectively a differentiation with respect to t and ^o- This phase velocity can be 
understood either as an angular velocity or a hnear tangential velocity along the swimmer surface (in radii per unit 
time). Equivalently, the wavelength A of the optimal surface deformation is defined as 



A = 2V0'i = 27r|14|, 



(68) 



as the wave period is equal to 27r. 

As the constraint on the maximum displacement Omax is stiffened (and the swimming efficiency is reduced), the 
magnitude of the wave velocity Vg and wavelength A are reduced (ip'i increas es). The modal decomposition is used 
to determine the phase velocity for all the optimal strokes presented in Sec. |III[ and the dependance of A (or Vg) 
with 6max is shown on Fig. 19 The maximum wave velocity is reached for the stroke of maximum efficiency (and 



displacement), for which the phase velocity magnitude is approximately equal to 0.67 (in radians) corresponding to a 
wavelength of about four radii (or two body lengths). For the stroke with smallest amplitude shown in Fig. 18 the 



phase velocity magnitude is reduced to about 0.27 corresponding to a wavelength of about 1.5 radii. Note that this 
reduction of the phase velocity and wavelengths with the stiffness of the constraint on the maximum displacement 
was also apparent by comparison of the slope of the recovery shock in Figs. |4]and[8j 

The phase velocity of the wave is negative in 9 (or positive in x) corresponding to a wave always traveling toward 
the north pole, in the same direction as the swimming velocity, and thus in the opposite direction to the effective 
stroke. Such metachronal waves are known as antiplectic [^I14j. This result is consistent with the Taylor's swimming 
sheet model for which the swimming velocity and the wave velocity have same direction for tangential deformations 
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FIG. 19: (color online) Variation of the wavelength A (or phase velocity IVe] = \/2ii) measured along the surface of the spherical 
swimmer, with the maximum angular displacement, Omax- The blue star corresponds to the unconstrained optimal swimming 
stroke of Sec. |II1 



[31 1341 I41j , whereas normal displacements lead to wave and swimming velocities of opposite directions (symplectic 
metachronal waves) [21 HZ] ■ 

Note finally that the waves quantified above are defined in terms of 9^ and t. The phase velocity is thus measured 
using the Lagrangian label rather than the actual instantaneous position of the cilia. Because the amplitude of beating 
is not small, this introduces a small difference between the quantitative measure of the phase velocity we propose and 
the phase velocity visually assessed, for example, from Fig. |4] A formal wave analysis in Eulerian coordinates leads 
to insignificant changes to the conclusions of the Lagrangian analysis. 



3. Mode amplitude 



and 14 
the 



The structure of the dominant mode presented on Fig. [18] also characterizes the amplitude of the tangential motion 
through the spatial and temporal amplitude functions ai(0o) and 6i(i). The temporal amplitude was normalized to 
unit mean value so that the spatial amplitude provides the physical value of the amplitude of displacement. We see 
that the departure from unity of hi is small confirming the structure of the dominant mode as a traveling wave with 
spatial variations of its amplitude. 

The spatial amplitude ai is observed to match very well the function Q{9o) obtained for each value of 6*0 as the 
amplitude of motion of the particular Lagrangian point. For the unconstrained stroke it is maximal at the equator 
As the constraint on maximum displacement is stiffened, the spatial amplitude is reduced (Omax in Figs, 
corresponds to the peak of this curve) and its profile is fiattened. As strokes with smaller Gmax are considered 
displacement amplitude of individual surface elements display smaller disparities between equatorial points and the 
rest of the surface. Due to the spherical symmetry, the shape function ai must go to zero near the poles and can 
therefore not be constant over the whole surface of the swimmer. However, we notice that for the most constrained 
case plotted in Fig. 18 the spatial amplitude is close to 8max for 30° < 9q < 150°. This result is particularly interesting 
as the strokes of ciliated microorganisms correspond to the lower end of the 0max-range considered here, where the 
shape function is the flattest. In such cases, the surface elements display, apart from the poles, a constant amplitude 
along the swimmer. This is consistent with the physical intuition that individual cilia do not "know" whether they are 
located near the pole or the equator and therefore, the beating properties (other than the phase) should be expected 
to be constant over the surface of the organism. 
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FIG. 20: (color online) Ratio of the swimming velocity to the wave speed, U/Vg, as a function of the maximum tangential 
displacement multiplied by the optimal wave number. As in other figures, results are binned and error bars are given both 
horizontally and vertically. The blue star corresponds to the unconstrained optimal swimming stroke of Sec. Ill] The two solid 
lines correspond respectively to U /Ve = (fe-Dmax)^/4 and U /Vg — (fc_Dmax)^/2 [3]. 



4- Scaling of the swimming-to-wave velocity ratio 



In the limit of small amplitude displacements, the swimming sheet model first introduced by Taylor shows that the 
swimming velocity (scaled by the wave propagation speed) is a quadratic function of the displacement amplitude [5ll27j. 
This scaling also applies to surface motions of larger amplitudes, and for a large variety of ciliated microorganisms, 
it has been established that U /c ~ a{kl)'^ with a ^ 0.25 ~ 0.5 and c is the wave velocity, I the cilia length and k the 
wave number (see Fig. 20 in Ref. [5). We plot in Fig. 20 the variations of CZ/IVgl with A:Z?niax where D^^ax = ©max 
is the maximum linear displacement amplitude and k — l/\Vg\ is the wavenumber of the optimal strokes (the stroke 
period and swimmer's radius have been normalized to 2tt and 1 respectively). We observe the same quadratic scaling 
of U with fcZ?,iiax, with a numerical prefactor which lies in the same range. 

A few comments must however be made. First, the results of Fig. 20 in Ref. [^ use the length I of an individual 
cilium rather than the maximum displacement. For most ciliated organisms, a good estimate of the cilia length is 
obtained as / ^ 2Dniax [!]• In addition, most ciliated organisms reviewed in Ref. [1] are not spherical but actually 
elongated in the swimming direction (such as the one shown in Fig. II]) . For such an organism, the wave velocity along 
the surface is very close to its component along the swimming direction, while these quantities differ by a factor 7r/2 
on a sphere. The factor U/Vx is therefore increased if the projection V^ of |Ve| along e^. is considered. As a result of 
these rescalings, one would find that U/c « 0.1 - 0.2 k^P, which is close but not exactly similar to the U/c « 0.25 - 
0.5 fc^/^ obtained in Ref. |4]. Additionally, we observe that the optimal strokes correspond to waves traveling faster 
than the swimming velocity, while many microorganisms show much slower wave speeds (resulting in [//jVel > 1). 
This limitation was already mentioned in Ref. |14) and points to the limits of the envelope model. 



VI. DISCUSSION 



In this paper we have considered a spherical envelope model (so-called squirmer) to investigate energetics in cilia 
dynamics and locomotion. Allowing only tangential but time-periodic deformations, we have used an optimization 
method based on a variational approach to derive computationally the stroke leading to the largest swimming efficiency. 
The optimal stroke was shown to display weak Lagrangian asymmetry, but strong Eulerian asymmetry, indicative of 
symmetry-breaking at the whole-organism level, but not at the level of individual cilia. We then added a constraint 
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in the optimization approach to penaUze large-amphtude deformation of the surface, together with a numerical 
ansatz, and derived a complete optimization diagram where all values of the swimming efficiency between and 50% 
(mathematical upper bound) could be reached. We were also able to construct a swimmer which is 50% efficient, 
although it is mathematically singular. The deformation kinematics of the optimal strokes were always found to be 
wave-like, which we analyzed using a formal modal decomposition approach, and is reminiscent of metachronal waves 
in cilia arrays. 

Our optimization procedure always leads to swimming strokes displaying antiplectic waves. This is, in fact, an 
intrinsic feature (and admittedly, a limitation) of the squirmer envelope model. In a real ciliated organism, the 
effective and recovery strokes correspond to very different shapes of the individual cilium [21 Ej . This asymmetry 
defines the effective and recovery strokes, independently from the coordination of neighboring cilia. In the present 
envelope model, such an asymmetry is not represented, as only the tangential displacement of cilia tips is prescribed, 
and thus the effective or recovery nature of the stroke is determined by the collective behavior. 

Previous calculations in the limit of small deformations and harmonic waves for the swimming sheet model ^ ^J 
or the squirmer model [34] have shown that the swimming velocity is always oriented in the same direction as the 
propagating wave. A physical argument leading to the same conclusion in the squirmer geometry can be proposed as 
follows. In the envelope model with tangential displacements, the half-stroke in the direction of the wave propagation 



corresponds to a compression of the surface. The definition of the average swimming velocity from Eq. (22) as a 
geometrically weighted-average of the Eulerian surface velocity identifies the effective stroke (of opposite direction to 
the swimming velocity) as the stroke where the surface is stretched because it occupies a greater part of the velocity 



map in the Eulerian (0,i)-plane (see Fig. 17 1). The wave and the effective stroke therefore have opposite directions, 
and the metachronal wave is antiplectic. 



This can also be seen mathematically. For the general swimming strokes in Eq. ( 64 1 , we can obtain using Eq. ( 22 ) 
that 

iU)^ve = -^f f {a{9Q)g'it-eo/vg)sm[eo+a{9o)g{t-eo/ve)]ydtdeo<0. (69) 



As a result, a deformation wave traveling from the south to the north pole {vg < 0) imposes a northward swimming 
velocity in the same direction as the wave. Although we do not prove this result here for arbitrary periodic functions 
i}{9o,t), it appears to be a general result. The simplified assumptions of the envelope model used here therefore 
prevent us to conclude on the relative efficiency of both types of metachronal waves. 

Our work was but a first attempt at an optimization approach to cilia dynamics, and could be extended in a 
variety of ways. First, and foremost, one should allow deformations normal to the swimmer surface to take place. 
In that case the swimmer would display time-periodic shape changes, and the spherical-harmonics framework will 
no longer be applicable. A full numerical approach would thus have to be implemented to derive the instantaneous 
swimming speed and energetics. Another possible extension would allow for non-axisymmetric deformation to occur, 
in which case the swimming kinematics would also be including a rotation. The dissipation arguments presented in 
Ref. [34| suggest that any such rotation would be sub-optimal, but it would be instructive to obtain that result from an 
optimization approach (see also [351132]). The straightforward extension of our results to other shapes with tangential 
deformations, in particular prolate spheroid-like, would also be biologically relevant. Finally, one should consider 
other biologically-relevant transport quantities, such as the flux of nutrients transported by the swimming-induced 
flow. 
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